home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / lisp / eulisp / mpfeel.lha / MPFeel / Plurals / Random / p_random.m < prev    next >
Text File  |  1992-04-27  |  2KB  |  72 lines
  1. /******************************************************************************
  2.  
  3.    A portable plural random number generator based on three linear
  4.    congruential generators.  -- Numerical Recipes in C pp. 210
  5.  
  6.    Copyright 1991,1992 Lee Iverson
  7.  
  8.    Author:
  9.      Lee Iverson <leei@mcrcim.mcgill.edu>,
  10.      McGill Research Centre for Intelligent Machines (McRCIM)
  11.  
  12.    See "copyright.h" for complete copyright information.
  13.  
  14. ******************************************************************************/
  15.  
  16. static char RCSid[] =
  17.   "$Id: p_random.m,v 1.2 92/01/23 17:34:05 leei Exp $";
  18.  
  19. #include <values.h>
  20. #include <math.h>
  21.  
  22. #include "p_randomI.h"
  23.  
  24. plural double
  25. p_frandom _IMPL_VOID()
  26. {
  27.     register plural double tmp;
  28.     register plural int j;
  29.  
  30.     if ( !rand_init__ ) init_p_random(0);
  31.  
  32.     ix1__ = (IA1*ix1__ + IC1) % M1;
  33.     ix2__ = (IA2*ix2__ + IC2) % M2;
  34.     ix3__ = (IA3*ix3__ + IC3) % M3;
  35.     j = 1 + ((97*ix3__)/M3);
  36.     tmp = r__[j];
  37.     r__[j] = (ix1__+ix2__*RM2)*RM1;
  38.     return tmp;
  39. }
  40.  
  41. plural double
  42. p_uniform_random _IMPL((lo, hi),
  43.                register plural double lo _AND
  44.                register plural double hi)
  45. {
  46.     return lo + p_frandom()*(hi-lo);
  47. }
  48.  
  49. plural double
  50. p_normal_random _IMPL((stdev), plural double stdev)
  51. {
  52.     static int iset = 0;
  53.     static plural double gset;
  54.     plural double fac, r, v1, v2;
  55.     plural double p_log(), p_sqrt();
  56.  
  57.     if ( !iset ) {
  58.     do {
  59.         v1 = p_uniform_random((plural double)-1.0,(plural double)1.0);
  60.         v2 = p_uniform_random((plural double)-1.0,(plural double)1.0);
  61.         r = v1*v1 + v2*v2;
  62.     } while ( r >= 1.0 );
  63.     fac = p_sqrt(-2.0*p_log(r)/r);
  64.     gset = v1*fac; ++iset;
  65.     return v2*fac*stdev;
  66.     } else {
  67.     --iset;
  68.     return gset*stdev;
  69.     }
  70. }
  71.  
  72.